/*unit.cpp*/
#pragma once
#include "unit.h"

unit::unit(node i, node j, node m) {
	list[0] = i.index();
	list[1] = j.index();
	list[2] = m.index();
	N_matrix[0] = j.x() * m.y() - j.y() * m.x();
	N_matrix[1] = j.y() - m.y();
	N_matrix[2] = m.x() - j.x();
	N_matrix[3] = m.x() * i.y() - m.y() * i.x();
	N_matrix[4] = m.y() - i.y();
	N_matrix[5] = i.x() - m.x();
	N_matrix[6] = i.x() * j.y() - i.y() * j.x();
	N_matrix[7] = i.y() - j.y();
	N_matrix[8] = j.x() - i.x();
	area = (N_matrix[0] + N_matrix[3] + N_matrix[6]) / 2;
}
//获取广义域序号在单元内的序号
int unit::find_index(int index) {
	for (int i = 0; i < 3; i++) {
		if (index == list[i]) return i;
	}
	return -1;
}
//狭义坐标形函数某一系数值
double unit::N_node_in(int inum, int flag) {
	return N_matrix[3 * inum + flag];
}
//单元内广义域序号对应节点形函数某一系数值(0=a,1=b,2=c)
double unit::N_node(int index, int flag) {
	int i = find_index(index);
	if (i == -1) return -1;
	return N_node_in(i, flag);
}
//单元内局部编号对应节点形函数在某点的值
double unit::N_point(int index, double x, double y) {
	int i = index;
	double N_xy = N_node_in(i, 0) + N_node_in(i, 1) * x + N_node_in(i, 2) * y;
	return N_xy / (2 * area);
}
//应变矩阵, 局部编号: f_div自由度(1=u,2=v); c_flag应变型(1=xx,2=yy,3=xy)
double unit::B_node(int index /*0,1,2*/, int f_div, int c_flag) {
	int i = index;
	double b = 0;
	if (f_div == 1) {
		if (c_flag == 1) b = N_node_in(i, f_div);
		else if (c_flag == 3) b = N_node_in(i, f_div + 1);
		else return 0;
	}
	else {
		if (c_flag == 2) b = N_node_in(i, f_div);
		else if (c_flag == 3) b = N_node_in(i, f_div - 1);
		else return 0;
	}
	return b / (2 * area);
}
//应力矩阵, 局部编号: f_div自由度(1=u,2=v); c_flag应力型(1=xx,2=yy,3=xy)
double unit::S_node(int index, int f_div, int c_flag, double E, double nu) {
	int i = index;
	if (c_flag == 3) f_div = f_div % 2 + 1;
	else if (f_div == 2) c_flag = c_flag % 2 + 1;
	double s = N_node_in(i, f_div);
	switch (c_flag)
	{
	case 1:
		s = s;
		break;
	case 2:
		s *= nu;
		break;
	case 3:
		s *= (1 - nu) / 2;
		break;
	}
	s = E * area * s / (2 * (1 - nu * nu));
	return s;
}
//刚度矩阵, f1,f2为自由度编号
double unit::K_node(int index1, int f1, int index2, int f2, double E, double nu, double t) {
	int i1 = find_index(index1);
	int i2 = find_index(index2);
	if (i1 == -1 || i2 == -1) return 0;
	double k = 0;
	if (f1 == 1) {
		if (f2 == 1) k = N_node_in(i1, 1) * N_node_in(i2, 1) + 0.5 * (1 - nu) * N_node_in(i1, 2) * N_node_in(i2, 2);
		else k = nu * N_node_in(i1, 1) * N_node_in(i2, 2) + 0.5 * (1 - nu) * N_node_in(i1, 2) * N_node_in(i2, 1);
	}
	else {
		if (f2 == 1) k = nu * N_node_in(i1, 2) * N_node_in(i2, 1) + 0.5 * (1 - nu) * N_node_in(i1, 1) * N_node_in(i2, 2);
		else k = N_node_in(i1, 2) * N_node_in(i2, 2) + 0.5 * (1 - nu) * N_node_in(i1, 1) * N_node_in(i2, 1);
	}
	k = E * t * 0.25 * k / ((1 - nu * nu) * area);
	return k;
}